Intro Questions

• Describe geometric anisotropy. Discuss this in the context of a data set.

• For Today:
  • Finish geometric anisotropy
  • Spatial Generalized Linear Models

Geometric Anisotropy Model

• Let \(\boldsymbol{Y}(\boldsymbol{s}) = \mu(\boldsymbol{s}) + w(\boldsymbol{s}) + \epsilon(\boldsymbol{s})\), thus \(\boldsymbol{Y}(\boldsymbol{s}) \sim N(\mu(\boldsymbol{s}), \Sigma(\tau^2, \sigma^2, \phi, B))\), where \(B = L^T L\).
• The covariance matrix is defined as \(\Sigma(\tau^2, \sigma^2, \phi, B)) = \tau^2 I + \sigma^2 H((\boldsymbol{h}^T B \boldsymbol{h}^T)^{\frac{1}{2}}),\) where \(H((\boldsymbol{h}^T B \boldsymbol{h}^T)^{\frac{1}{2}})\) has entries of \(\rho((\boldsymbol{h_{ij}}^T B \boldsymbol{h_{ij}}^T)^{\frac{1}{2}}))\) with \(\rho()\) being a valid correation function, typically including \(\phi\) and \(\boldsymbol{h_{ij}} = \boldsymbol{s_i} - \boldsymbol{s_j}\).

Geometric Anisotropy Visual

• Consider four points positioned on a unit circle.

• How far apart are each set of points?

Geometric Anisotropy Exercise 1

Now consider a set of correlation functions. For each, calculate the correlation matrix and discuss the impact of \(B\) on the correlation. Furthermore, how does B change the geometry of the correlation?

1. \(\rho() = \exp(-\boldsymbol{h_{ij}}^T B \boldsymbol{h_{ij}}^T)^{\frac{1}{2}})),\) where \(B = \begin{pmatrix} 1 & 0 \\ 0 & 1 \\ \end{pmatrix}\)

2. \(\rho() = \exp(-\boldsymbol{h_{ij}}^T B \boldsymbol{h_{ij}}^T)^{\frac{1}{2}})),\) where \(B = \begin{pmatrix} 2 & 0 \\ 0 & 1 \\ \end{pmatrix}\)

3. \(\rho() = \exp(-\boldsymbol{h_{ij}}^T B \boldsymbol{h_{ij}}^T)^{\frac{1}{2}})),\) where \(B = \begin{pmatrix} 3 & 1 \\ 1 & 1 \\ \end{pmatrix}\)

Geometric Anisotropy: Solution 1

\(\rho() = \exp(-\boldsymbol{h_{ij}}^T I \boldsymbol{h_{ij}}^T)^{\frac{1}{2}}))\)

Implied Distance

Correlation

Geometric Anisotropy: Solution 2

\(\rho() = \exp(-\boldsymbol{h_{ij}}^T B \boldsymbol{h_{ij}}^T)^{\frac{1}{2}})),\) where \(B = \begin{pmatrix} 2 & 0 \\ 0 & 1 \\ \end{pmatrix}\)

Implied Distance

Correlation

Geometric Anisotrop: Solution 3

\(\rho() = \exp(-\boldsymbol{h_{ij}}^T B \boldsymbol{h_{ij}}^T)^{\frac{1}{2}})),\) where \(B = \begin{pmatrix} 3 & 1 \\ 1 & 1 \\ \end{pmatrix}\)

Implied Distance

Correlation

More Geometric Anisotropy

• The matrix \(B\) relates to the orientation of a transformed ellipse.
• The (effective) range for any angle \(\eta\) is determined by the equation \[\rho(r_\eta(\tilde{\boldsymbol{h}}_{\eta}^T B \tilde{\boldsymbol{h}}_{\eta}^T)^{\frac{1}{2}}) = .05,\] where \(\tilde{\boldsymbol{h}}_{\eta}\) is a unit vector in the direction \(\eta\).

Fitting Geometric Anisotropy Models

• Okay, so if we suspect that geometric anisotrophy is present, how do we fit the model? That is, what is necessary in estimating this model?
• In addition to \(\sigma^2\) and \(\tau^2\) we need to fit \(B\).
• What about \(\phi\)? What is \(\phi\) when \(B = \begin{pmatrix} 1 & 0 \\ 0 & 1 \\ \end{pmatrix}\)?

Priors for B

• While \(B\) is a matrix, it is just another unknown parameter.
• Hence, to fit a Bayesian model we need a prior distribution for \(B\).
• One option for the positive definite matrix is the Wishart distribution, which is a bit like a matrix-variate gamma distribution.

Data Viz Motivation

Ozone Exceedance in Colorado

Generalized Linear Model Notation

There are three components to a generalized linear model:

1. Sampling Distribution: such as Poisson or Binomial
2. Linear combination of predictors: \(\eta = X\beta\)
3. A link function to map the linear combination of predictors to the support of the sampling distribution.

Logistic Regression Overview

Write out the complete model specification for logistic regression.

• Assume \(Y_i\) is the binary response for the \(i^{th}\) observation, \[\begin{eqnarray*} Y_i &\sim& Bernoulli(\pi_i)\\ logit(\pi_i) &=& X_i \beta, \end{eqnarray*}\]
• where \(logit(\pi_i) = log \left(\frac{\pi_i}{1-\pi_i}\right)\)

glm()

Interpret this output

CO <- CO %>% mutate(north = as.numeric(Latitude > 38 ))
glm(Exceedance~north, family=binomial(link = 'logit'),data=CO) %>% summary

## 
## Call:
## glm(formula = Exceedance ~ north, family = binomial(link = "logit"), 
##     data = CO)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -2.313   0.378   0.378   0.378   1.177  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)  
## (Intercept) -1.695e-15  8.165e-01   0.000   1.0000  
## north        2.603e+00  1.097e+00   2.372   0.0177 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 28.708  on 34  degrees of freedom
## Residual deviance: 22.873  on 33  degrees of freedom
## AIC: 26.873
## 
## Number of Fisher Scoring iterations: 5

Logistic Regression in JAGS

JAGS specification for logistic regression

model{
  # likelihood
    for (i in 1:N){
        y[i] ~ dbern(p[i])
        p[i] <- 1 / (1 + exp(-eta[i]))
        eta[i] <- int + x[i] * north
    }
    # prior
    int ~ dnorm(0, 1E-12)
    north ~ dnorm(0, 1E-12)
}

Spatial Logistic Regression

• Assume \(Y(\boldsymbol{s_i})\) is the binary response for \(\boldsymbol{s_i}\), \[\begin{eqnarray*} Y(\boldsymbol{s_i})|\beta, w(\boldsymbol{s_i}) &\sim& Bernoulli(\pi(\boldsymbol{s_i}))\\ logit(\pi(\boldsymbol{s_i})) &=& X(\boldsymbol{s_i}) \beta + w(\boldsymbol{s_i}), \\ \end{eqnarray*}\]
• where \(\boldsymbol{W} \sim N(\boldsymbol{0},\sigma^2 H(\phi))\)

Spatial Logistic Regression in JAGS

logistic.spatial <- "model {
  for (i in 1:N) {
      Y[i] ~ dbern(p[i])
      logit(p[i]) <- mu[i] + w[i]
      mu[i] <- beta[1]+beta[2]*x.north[i]
    muW[i] <- 0
  }
  # process
  w[1:N] ~ dmnorm(muW[],Omega[,])

  # priors
  beta[1] ~ dnorm(0.0,1E-8)
  beta[2] ~ dnorm(0.0,1E-8)
  spat.prec ~ dgamma(0.001, 0.001)
  sigmasq <- 1/spat.prec
  phi ~ dunif(1.5,2.5)  
    
  #build omega
  for (i in 1:N){
    for (j in 1:N){
      H[i,j] <- (1/ spat.prec) * exp(-phi * d[i,j])
    }
  }
  Omega[1:N,1:N] <- inverse(H[1:N,1:N])
}"

spGLM()

m.1 <- spGLM(y~1, 
             family="binomial", 
             coords=coords, 
             weights=weights, 
             starting=list("beta"=beta.starting, "phi"=0.06,"sigma.sq"=1, "w"=0),
             tuning=list("beta"=beta.tuning, "phi"=0.5, "sigma.sq"=0.5, "w"=0.5),
             priors=list("beta.Normal"=list(0,10), "phi.Unif"=c(0.03, 0.3), "sigma.sq.IG"=c(2, 1)),
             amcmc=list("n.batch"=n.batch, "batch.length"=batch.length, "accept.rate"=0.43),
             cov.model="exponential", 
             verbose=TRUE, 
             n.report=10)

Motivation

Poisson Regression Overview

Write out the complete model specification for Poisson regression.

• Assume \(Y_i\) is the count response for the \(i^{th}\) observation, \[\begin{eqnarray*} Y_i &\sim& Poisson(\lambda_i)\\ \log(\lambda_i) &=& X_i \beta, \end{eqnarray*}\]
• thus \(\exp(X_i \beta) \geq 0\)
• next write out a Poisson regression model with spatial random effects

Spatial Poisson Regression

• Assume \(Y_i\) is the count response for \(\boldsymbol{s_i}\), \[\begin{eqnarray*} Y(\boldsymbol{s_i}) &\sim& Poisson(\lambda(\boldsymbol{s_i}))\\ \log(\lambda(\boldsymbol{s_i})) &=& X(\boldsymbol{s_i}) \beta + w(\boldsymbol{s_i}), \end{eqnarray*}\]
• where \(\boldsymbol{W} \sim N(\boldsymbol{0},\sigma^2 H(\phi))\)

Fitting Spatial Poisson Regression Models

1. spGLM() accepts binomial and Poisson data,
2. JAGS can also be used to fit spatial GLMS models, or
3. Try writing a sampler from scratch.